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ABSTRACT 

Current theories on planetary formation establish that giant planet formation 
should be contextual to their quick migration towards the central star due to the 
protoplanets-disc interactions on a timescale of the order of 10 5 years, for objects 
of nearly 10 terrestrial masses. Such a timescale should be smaller by an order of 
magnitude than that of gas accretion onto the protoplanet during the hierarchical 
growing-up of protoplanets by collisions with other minor objects. These arguments 
have recently been analysed using N-body and/or fluid-dynamics codes or a mixing of 
them. In this work, inviscid 2D simulations are performed, using the SPH method, to 
study the migration of one protoplanet, to evaluate the effectiveness of the accretion 
disc in the protoplanet dragging towards the central star, as a function of the mass of 
the planet itself, of disc tangential kinematics. To this purpose, the SPH scheme is con- 
sidered suitable to study the roles of turbulence, kinematic and boundary conditions, 
due to its intrinsic advective turbulence, especially in 2D and in 3D codes. Simulations 
are performed both in disc sub-Keplerian and in Keplerian kinematic conditions as 
a parameter study of protoplanetary migration if moderate and consistent deviations 
from Keplerian Kinematics occur. Our results show migration times of a few orbital 
periods for Earth-like planets in sub-Keplerian conditions, while for Jupiter-like plan- 
ets estimates give that about 10 4 orbital periods are needed to half the orbital size. 
Timescales of planet migration are strongly dependent on the relative position of the 
planet with respect to the shock region near the centrifugal barrier of the disc flow. 

Key words: planetary systems: formation - planetary systems: protoplanetary discs 



1 INTRODUCTION 

Currently, nearly 300 extra-solar planets have been detected 
(http://www.exoplanets.eu/), and most of them are mas- 
sive (Jupiter-like) planets with small (fractions of AU ) semi - 
major axes ( see lPerrvmar] <|2000f > and lUdrv fc Santos! ([2007) 
for a review on detection methods and on the main prop- 
erties of exo-planets) . Theoretical models for the formation 
of planetary systems are then mostly based on the crucial 
role of an accretion disc around a forming young star, where 
the disc provides both the material for the forming plan- 
ets and the interactions responsible for the migration of the 
protoplanets towards the central star. The need for proto- 
planetary migration comes from the difficulty of an "in situ" 
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formation in such close to the star positions, so that the 
detailed understanding of migration mechanisms appears 
mandatory in order to build self-consistent models for plane- 
tary systems formation. Some analytic studies on the matter 
llGoldreich fc Tremaine ! ll979l . ll98ol : iLin fc Papaloizoulll984l . 
1986al lbl: IWardlll997f) suggest that gravitational torques due 
to resonant non-axisymmetric structures could act to trans- 
fer angular momentum between the protoplanets and the 
disc ( for a discussion o f other alternative migration mod- 
els see iDel Popolo et all (120051 ) and references therein) . The 
gravitational disturbance exerted by the protoplanet on the 
disc produces spiral density waves which break the axial 
symmetry of the disc structure. Perturbation theory allows 
an estimate of the net torque that this perturbed disc ex- 
erts back onto the planet at "Lindblad resonance sites" 
l|Goldreich fc Tremaindll98d : IWardlll997l ). This disc-planet 
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interaction induces also a repulsion of material on either 
side of the protoplanet's orbit, with a possible formation 
of a density gap, depending on the interaction strength. 
Two main migratio n types are then distinguished: type I 
|Tanaka et alj 12002). when an Earth-like planet, not capa- 
ble of opening a density gap in the disc is embedded in 
the disc structure. In this case the migration rate should 
be proportional to the protoplanet mass; type II when a 
massive protoplanet, capable of opening and maintaining 
a density gap, remains trapped in this gap. In this case 
migration is due to the viscosity dominated evolution of 
the disc. In both cases it is possible that the lifetime of 
migration is shorter than the lifetimes of planet forma- 
tion or disc dissolution, causing the collapse of the planet 
toward the star. Accurate numerical sim ulations are nec- 
essar y to better understand the matter l|Papaloizou et all 
2006). The main aim of this work is to study the role 
of initial and boundary conditions on planetary migration 
as a consequence of planet-disc dynamic interaction dur- 
ing the initial formation phases of stellar systems. Several 
papers devoted to this theme, according to various fluid- 
dynamics schemes jArtvmowi cz 2004; D'Angelo et al.l 20021. 
2001 120061 ; lKley||2000l ; iPapaloizou et al.ll2006l ; ISchafer et al l 
2004 ). showed that this mechanism should be very effec- 
tive for the protoplanet dragging toward the central proto- 
star in characteristic time-scales of the order of 10 5 years. 
These studies assume that the accretion disc is basically 
Keplerian and pressure forces are often neglected among 
the protostar-disc interactions. In Keplerian models, the 
Keplerian distribution of disc g as velocity is usually as- 
sumed as a reasonable hyp othesis l|Cresswell fc Nelson! 
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ypothcsis (Orcsswcll & INclso 
jj iKlevI |2000| ; ISchafer et al 
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Terauem et aiT feoOO) . However, as other authors have al- 



ready done for a wider context du r ing the accre tion disc for- 
mation phase (|Sollins fc Holl2005l ; IUlrichlll976l ) or for small 
sub-Kepler ian correction due to the contrib ution of pressure 
gradients ( Paa rdekooper fc Mellemal 120041 ) . we find useful 
to explore here the role of sub-Keplerian flows as an open 
possibility for pre-main sequence accretion, as well as their 
role in the radial transport of solid bodies. Therefore sub- 
Keplerian flows in the formation of protoplanetary accretion 
discs are here analysed as a parametric study (see section 
2.21 for details on exact initial conditions). 
In this paper non viscous Cartesian 2D SPH sub-Keplerian 
disc models are given, with planetary migration times both 
for Jovian and for Terrestrial planets. Models of Jovian or 
Terrestrial planets, interacting with an initial Keplerian ac- 
cretion disc, are also produced for a strict comparison. 
The coming out of shock fronts and complex structures, both 
steady and progressive in sub-Ke plerian turbulent accre- 
tion fl o ws in AGN, wa s descr ib ed inlChakrabarti fc Moltenil 
|l993t ). iMolteni et ail |l994h . ILanzafame et all ll 19981), as 
well a s pulsed or periodic outflows in ILanzafame et al.l 
(2008) as a consequence of coUisional interactions in the disc 
surrounding the central accreting object in axial symmetry. 
In such environments (AGN or equivalently protostellar sys- 
tems) , the potential importance of the centrifugal barrier in 
developing outflow structures in turbulent coUisional fluid- 
dynamics, was successfully described in those papers. This 
implies that the possibility of a slower dragging toward the 
central star cannot be excluded in principle. Therefore, as a 
secondary aim, the purpose of the paper is the methodical 



search for those initial and boundary conditions determin- 
ing results showing how Jovian and Terrestrial protoplanets 
could stay in orbits surrounding the central stars at dis- 
tances of the order of that of Venus or Mercury. 

In section [2] we discuss the initial and boundary con- 
ditions characterising our models. In sec. Owe describe our 
results, and we discuss them in section [4] 



2 MODELS FEATURES AND BOUNDARY 
CONDITIONS 

2.1 Models 

The simulations presented in this work are based on a two- 
dimensional model which makes use of the Smo othed Par- 
ticle Hydrodynamics (SPH) numerical method l|Monaghanl 
1992). The choice of using the SPH method is funded on the 
grounds of its Lagrangian nature and its ability to easily 
tackle hydrodynamic problems with free surfaces. The abil- 
ity of many different numerical schemes, both grid-based and 
Lagrangia n, to correctly descr i be pla net-gas interactions was 
tested by Ide Val-Borro et all l|200cJ) . showing that SPH is 
not significantly less performant than other methods. Given 
the lack of a wide range of case studies on the matter in SPH, 
we believe that our work could give an interesting contri- 
bution to improve the basic knowledge on disc-protoplanet 
interaction. 

The code does not contain a specific radiation treatment, 
but the used equation of state: P = (7 — l)pit, where u is 
the thermal energy per unit mass, includes the polytropic 
index 7, which can be adjusted to values lower than 5/3 if 
radiation, partial molecular dissociation or partial ionisation 
effects are present, giving the overall effect of a higher com- 
pressibility to the gas. Current calculations are performed 
with the simple choice of 7 = 5/3, but we intend to further 
investigate this matter in a future work. 
The model is built in Cartesian coordinates with dimension- 
less quantities. Reference physical units are: 

- the initial stellar mass Mo; 

- the stellar radius 7?o; 

- the Keplerian period of an orbit of radius Ro around a 



star of mass Mo, To 



-R 



3/2 



With these choices, the following units can be used to trans- 
form our numerical results in physical quantities: 



GMg 

R 



velocity: vo = 

specific angular momentum: jo 
energy per unit mass: 



Mo 



surface density: £0 = -gf 



= ^VGMoRo 

GMg 



Ro 



The initial conditions used to build the accretion disc are ax- 
ially symmetric. The accretion disc is in fact created through 
the injection of particles at point-like positions (injectors) 
along circles concentric with the central star. Particles have 
initial tangential velocities set by choosing a value for the 
angular momentum per unit mass j. At some point, the ac- 
cretion disc reaches a "nearly" steady state population of 
particles (time averaged since some variation of the number 
of particles and of the density distributions are actually ob- 
served during our simulations, due to the highly turbulent 
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nature of our calculated flows) with 3x 10 5 — 5x 10 5 particles, 
depending on the model, when using a smoothing length 
(spatial resolution parameter of the SPH method) h = 0.3. 
At this stage, the planet is inserted and the model is evolved. 
We included the following interactions in our models: 

(i) gravitational interaction between gas particles and 
central star; 

(ii) gravitational interaction between planet and star; 

(iii) gravitational interaction between particles and 
planet; 

(iv) gas pressure between neighbouring gas particles; 

(v) artificial viscosity (pressure contribution) between gas 
particles; 

(vi) gas particle capture by the planet (conservation of 
momentum is used to correct the planet speed before re- 
moving the particle) when gas particles approach the planet 
area (planet-particle distance less than 2h). 

Additionally, in some of our models, we tried to add the 
above missing pressure interaction between the disc gas and 
the protoplanet. With the above specified conditions, our 
models do not have the necessary spatial resolution to treat 
the details of the interaction between the gas and the proto- 
planet surface. The simplest way to add pressure interaction 
in this case, is to add a further "ad hoc" SPH gas particle 
at the protoplanet position, giving to it the planet's veloc- 
ity and the mass of the other surrounding gas particles. The 
density is then calculated according to the SPH method and 
the pressure force acting onto this particle is directly trans- 
ferred to the planet itself. Though this approach, which is 
hereafter indicated as adding a "pseudo-atmosphere" to the 
protoplanet, may appear simplistic, it allows the introduc- 
tion of the, otherwise totally missing, pressure interaction 
acting onto the planet. 

No self-gravitation is included in the gas for computational 
reasons, and this is not believed to negatively affect our com- 
putations since the disc total mass is 10 -11 x the number of 
particles, which for ~ f 5 particles means fO" 6 (in units of 
the central star mass Mo), so that the dominating gravita- 
tional force is still that of the central star. Reported val- 



ues for disc masses are in the range fO 



ftT 1 M on 



the grounds of observational constrain ts from T-Tauri stars 
l|Hartmannl2000l : lTerquem et al.ll2000l ). but this value is usu- 
ally distributed over a disc with a radius of 100 — 1000 AU. 
In order to compare this to the structure of our simulated 
disc, given the dimensionless nature of our simulation, wc 
need to choose some physical reference unit s. If 1 Mp an d 
2-3 Pi,© are taken for a typical T-Tauri star (Bertout 1989), 
the disc has a dimension of about ~ 2 AU. Therefore our 
simulations cover only the central part of a typical T-Tauri 
accretion disc, and the mass included in the simulation is in 
agreement with observational constraints. 

Most accretion disc models in lit erature include vis- 
cosity , according to the disc model of IShakura fc Sunvaevl 
1 19731 ). Even today, it is hard to find a self-consistent way 
to give a correct value for the ass viscosity parameter. 
We choose instead to not include physical viscous terms in 
the equations of our models, because we are here focused 
on the protoplanet migration mechanism based on gravita- 
tional and pressure forces and momentum transfer through 
gas particle captures, rather than to viscosity driven migra- 
tion. Our models our non viscous, meaning that ar based on 




60 
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Figure 1. Azimuthally averaged density profiles of the sub- 
Kcplcrian accretion discs, just before planet insertion, for j = 
18, 36, 54. 



the Euler equations. The absence of physical viscosity en- 
hances the turbulent nature of our simulation, allowing us 
to correctly describe some effects due to turbulence like the 
development of non-axisymmetric structures in the disc (in- 
duced also by the presence of the protoplanet). Turbulence 
might also produce stronger shock fronts when infalling gas 
particles hit the centrifugal barrier. 

There are plans to include physical viscosity in a future 
work, in order to better explore the viscosity driven migra- 
tions (like in type II migration models) and viscosity contri- 
bution in general. 



2.2 Initial conditions 

Three choices are made here for the angular momentum 
per unit mass j: 18, 36, 54, for particles injected at a radial 
distance of 130. All the three choices give a sub-Keplerian 
feature to the velocity distribution of the accretion disc 
(particles injected at a radial distance of 130 would have 
Keplerian velocity with j ~ 70), so that injected particles 
start decreasing their distance to the centre, until they "hit" 
the centrifugal barrier. These choices are considered in the 
the framework of a parametric study, to test the effect of a 
wide range of possible conditions, from "extreme" (j = 18) 
to "weak" (j — 54) sub-Keplerian. The mass of the particles 
is set to 10 -11 (in units of Mo), which means that the total 
mass of the disc can be about 10 -6 — 10 -5 . Calculations are 
performed with two masses for the planet: 10 -7 (Earth-like 
planet) and 10 -3 (Jupiter- like planet). The computational 
domain is circular, extended in the xy plane with a radius 
of 150 (the star is located in the "origin"), and the planet 
is located inside the disc plane, with the initial position at 
X p = 100 and Y p = 0, and an initial zero eccentricity. 

The azimuthally averaged density profiles of the sub- 
Keplerian accretion discs, just before planet insertion, are 
plotted in figure Q] for the three cases (j = 18,36,54). A 
Density peak is present near the centrifugal barrier in each 
case. For j = 18 the peak is close to the star, while for 
higher j values the density close to the star has values 
similar to those found at the disc outer edges. Density in 
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36b 
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Table 1. List of computed models, with their identifying labels 
(first column). The second column gives the specific angular mo- 
mentum of injected gas particles, the third column gives the ini- 
tial mass of the protoplanet and the fourth column indicates the 
eventual presence of the pseudo-atmosphere on the protoplanet. 



SPH is computed according to local particle mass distri- 
bution. At the disc outer edge, particles are produced and 
injected in the stellar gravitational potential well. Instead 
particles are in a free fall very close to the central star, 
where gravitational potential rules the whole gas dynamics, 
and they are "removed" from the simulation (the mass of 
the particle is added to the star mass) when they reach a 
radial position equal to the star radius. 
As far as the disc material accretion rate onto the central 
star, if 1 Mb and 2-3 R© are taken for a typical T-Tauri 
star <|Bertoutill989h . we obtain for the three cases values of 
2 x IO" 6 , 6 x IO" 7 , 3 x IO" 7 M Q /year with oscillations of 
about 30%, 80%, 100%. These values are in agreement with 
estimated accretion rates for T-Tauri stars jGullbring et al.l 
1998; lHartigan. Edwards fc Ghandourll 19951 ). 
In order to compare some of our results with more tradi- 
tional models, some computations (two with an Earth-like 
planet and two with a Jupiter-like planet) are built using 
an initial exactly Keplerian distribution of gas particle 
velocities, with an initial uniform density over the entire 
disc area. For simplicity, we will refer to these models 
as "Keplerian", even if the velocity distribution of the 
gas particles cannot remain exactly Keplerian during the 
simulation, due to the pressure forces within the gas and to 
the interactions between the planet and the accretion disc. 
These models will be used to test the eventual formation of 
a gap near the protoplanet lo cation, and to compare some 
classical views on the subject <|Wardlll997n with our results 
obtained with an inviscid model. 

A summary of the distinguishing features of our models 
is given in table Q] together with identifying labels. In the 
table, M p is the mass of the planet in units of the central 
star mass, and "pseudo-atmosphere" indicates the eventual 
presence of a gas SPH particle associated with the planet. 



3 RESULTS 

The turbulent nature of our disc models is evidenced by 
the high value of the Reynolds number, eval uated by using 
the a rtificial viscosity coefficient u„ um = c s h |Molteni et al.l 
1991), especially at the disc outer edge, where we have 
values of about 1200, 2400, 3600 for j = 18, 36, 54 respec- 
tively. Instead, in the innermost regions, where the disc gas 
temperature is higher, we still evaluate Re values of about 
300 + 1000. 

Figures [2]|9] show the evolution of the orbital parameters 
for all models. The semi-major axis a as a function of 
time is shown in the top rectangle of each figure; the 
angular momentum of the protoplanet per unit mass l z , 
the coordinate x and the eccentricity e are shown in the 
following rectangles, respectively. The coordinate x is 
plotted to show the orbital phase, which allow us to count 
protoplanet orbits. Each figure includes the results from 
two models, with (models with labels "b" and "d") and 
without (models "a" and "c") the pseudo-atmosphere gas 
particle for the planet. The main feature that emerges 
from the sub-Keplerian models is a downward migration 
of the protoplanet, the speed of which being strongly 
dependent on the planet mass. Earth-like planets migrate 
in a few orbits, especially in the j — 18 case, while for the 
Jupiter-like planets, an extrapolation of our results (which 
are limited to 20-30 orbits), suggests that a few thousands 
orbits are necessary to move the planets towards the central 
star. 

Of course, the above suggested extrapolation would require 
a linear behaviour of the semi-major axis (or the specific 
angular momentum of the planet) with respect to time. 
But the results obtained for the less massive planet suggest 
that the migration rate asymptotically decreases for all the 
three values of the sub-Keplerian angular momentum of the 
injected particles. All the sub-Keplerian models with the 
Jovian planet lack the "classical" gap opening around the 
planet. This is not believed to be a problem as the type II 
migration mechanism, with the planet blocked inside the 
induced density gap, was developed in the framework of a 
Keplerian gas distribution in the accretion disc and within 
a viscosity based fluid model. Therefore the migration 
mechanism suggested within our model cannot be directly 
associated to the classical type I - type II distinction. 
For the Keplerian models, fig. [S] shows that the Jupiter-like 
planet migration slowly proceeds outwards, despite the 
development of density gap soon after the introduction of 
the planet in the initially totally Keplerian gas distribution 
(fig. I16[) . while the Earth- like planet (fig. [9} migrates 
downward, with a slower rate when compared to that of 
sub-Keplerian discs. We have in fact 0.6% variation in a in 
about 30 orbits, while, with the same number of orbits, the 
sub-Keplerian discs show a ~ 90%; 80%; 50% decrease of 
the planet semi-major axis, for j z — 18; 36; 54 respectively. 
As far as the influence of the pseudo-atmosphere is con- 
cerned, the difference between the results obtained with 
and without this model structure is of the order of ~ 10~ 2 % 
with the Jovian protoplanet and of ~ 1 — 4% with the 
Terrestrial protoplanet. In most cases we have an initial 
faster migration rate with the pseudo-atmosphere, but the 
situation is inverted after a few protoplanet's orbits. This is 
clearer in models with a Terrestrial protoplanet for which 
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Figure 2. Orbital parameters evolution for models 18a and 18b. 
The semi-major axis a, the specific angular momentum l z a, the x 
coordinate of the protoplanet and the eccentricity e are reported, 
as a function of time. The dotted and continuum graphs are for 
the two reported data sets, with and without the protoplanet 
pseudo-atmosphere. 



migration rate is much faster in all models. A significant 
difference is registered between the 54c and 54d models, 
where a ~ 10— 15% difference is obtained in the semi-major 
axis a in advanced migration phases. The model without 
the pseudo-atmosphere (54c) suffers the fastest migration 
rate at some point. 

In order to draw a semi-analytic behaviour in the orbital 
parameters evolution of the planets in our models, best-fits 
have been computed for the specific angular momentum j 
(figures I10I13[) to have an estimate of the time derivatives 
and decrease lifetime of the angular momentum. Two kind 
of fits have been used: linear and exponential, according to 
the following schemes: 

j z — a * exp(b ■ t) + c * exp(d ■ t) + e (1) 

jz =at + b (2) 

Tables [5] and report the best-fit parameters obtained for 
the various models. 

In table[2]the last two columns give the inverse values of 
the parameters b and d, and they have been dimensioned in 
years by choosing 1 Mg and 2 Rq as reference units. Of the 
two time-scales estimated in this way, the most important 
is the one with the corresponding highest "multiplying 
factor for the exponential" (a or c). For Jupiter- like planets 
these estimates must be taken with some care, as the 
simulations did not go enough far in time to have a clear 
path of the migration process, so that migration timeline of 
these models could be different than what expected from 
extrapolations. 

In table [3] the 4th column gives a dimensionless migration 
time, evaluated as l z /l z , and the last column gives the same 



Figure 3. Orbital parameters evolution for models 18c and 18d. 
See fig. [2] caption for description. 



j= 36 
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Figure 4. Orbital parameters evolution for models 36a and 36b. 
See fig. [2] caption for description. 



time converted in years by choosing the same reference 
units as in tabled 

For models a and b (Jovian planet) a linear fit already gives 
a good description of the angular momentum evolution of 
the planet, but an exponential fit has also been performed 
because, on a longer time-scale (a few thousands orbits), 
the curve showing j as a function of time could resemble 
that of the c and d models. We are therefore able to analyse 
the different migration lifetimes through comparable pa- 
rameters (exponential lifetime or linear costant derivative). 



6 V. Costa, V. Pirronello, G. Belvedere, A. Del Popolo, D. Molteni, G. Lanzafame 



IIIU-LLL.! 




|-> 












\ (hi i ~\ rci v G i 


"1 / fi ( ~\ /"O OT'C 1 

J- / LL \ y l_> cLl a 1 


18a 


62.87 


-3.4 x 10~ 8 


-3.412 x 10 


-2 


1.059 x 10" 5 


2.2 x 10" 


2 


-2.66 x 10 4 


8.54 x 10 


18b 


-1.054 x 10~ 3 


7.134 x 10" 5 


62.84 




-3.994 x 10~ 8 


0.0 




1.27 x 10 


-2.27 x 10 4 


18c 


138.4 


-4.426 x 10" 4 


7.672 




3.873 x 10- 5 


0.0 




-2.04 


2.34 x 10 


18d 


133.1 


-4.547 x 10" 4 


10.12 




1.822 x 10~ 5 


0.0 




-1.99 


4.97 x 10 


36a 


62.84 


-3.02 x 10" s 


1.228 x 10" 


-4 


1.519 x 10~ 4 


-4.0 x 10 


-3 


-3.00 x 10 4 


5.96 


36b 


62.83 


-3.209 x 10~ 8 


1.538 x 10" 


-3 


8.522 x 10~ 5 


5.0 x 10- 


3 


-2.82 x 10 4 


1.06 x 10 


36c 


77.65 


-1.842 x 10~ 4 


7.86 




6.437 x 10~ 5 


0.0 




-4.91 


1.41 x 10 


36d 


63.92 


-2.191 x 10~ 4 


20.37 




7.021 x 10~ 6 


0.0 




-4.13 


1.29 x 10 2 


54a 


2.84 


-1.557 x 10" 8 


7.47 x 10- 


6 


2.087 x 10~ 4 


-4.0 x 10 


-3 


-5.81 x 10 4 


4.33 


54b 


0.1758 


-1.734 x 10~ 5 


62.66 




2.504 x 10~ s 


2.2 x 10- 


3 


-5.22 x 10 


3.61 x 10 4 


54c 




















54d 


56.29 


-4.93 x 10~ 5 


14.5 




3.032 x 10~ 5 


0.0 




-1.84 x 10 


2.98 x 10 


kcpLa 


62.83 


5.484 x KT 10 


6.264 x 10" 


-5 


8.953 x 10~ 5 


1.7 x 10" 


3 


1.65 x 10 6 


1.01 x 10 


kepLb 


62.83 


4.825 x 10" 10 


1.375 x 10" 


-4 
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1.87 x 10 6 


1.36 x 10 


kcpLc 


65.74 


-6.256 x 10~ 7 


-2.903 




-1.407 x 10~ 5 


0.0 




-1.45 x 10 3 


-6.43 x 10 


kcpLd 


-1.982 


1.227 x 10- 5 


64.83 




3.28 x 10~ 7 


0.0 




7.37 x 10 


2.76 x 10 3 



Table 2. Exponential best-fit parameters for the specific angular momentum time evolution. A rough estimate of migration times is 
given by the inverse of the two parameters b and d, dimensioned values in years are given in 7th and 8th columns, calculated under the 
hypothesis of 1 Mq star with a radius of ~ 2 Rq . 
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time 



Figure 5. Orbital parameters evolution for models 36c and 36d. 
See fig. [2] caption for description. 
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Figure 6. Orbital parameters evolution for models 54a and 54b. 
See fig. [2] caption for description. 



For some of c and d models, the behaviour is initially linear 
and it becomes more exponential once the protoplanet 
reaches a position closer to the star. In these cases, the 
linear fit has been done for a limited time interval during 
which the evolution does not show a significant migration 
rate variation. A peculiar case is the 54c model where an 
initial linear behaviour with negative derivative is followed 
by a nearly complete migration stop. The exponential fit 
has not been evaluated for this model. 

Propagation of waves coming from the central zone of 
the disc has been observed, at least for the j = 18 and j — 36 
models. Some snapshots showing this phenomenon are plot- 
ted in figures [14] and [15] The same is not clearly observed 



in the other models, the cause being probably the high an- 
gular momentum of particles. A strong sub-Keplerian con- 
dition for injected particles makes them rapidly fall towards 
the centrifugal barrier, and the subsequent bouncing-back 
could be at the origin of these waves. Details on this prop- 
agation of shock waves i n systems with a centra l accreting 
object was described by lLanzafame et all (|2008h and refer- 
ences therein. A proper combination of sub-Keplerian gas 
flow and viscosity parameter values can create a complex 
shock structure with nearly periodical radial motions of the 
shock front. Periods of quiet piling up of particles near and 
below the centrifugal barrier are followed by the develop- 
ment of a shock front propagating outwards. In low or no 
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j= 54 




j= kepi 




Figure 7. Orbital parameters evolution for models 54c and 54d. 
See fig. [2] caption for description. 



Figure 9. Orbital parameters evolution for the Keplerian models 
and the Earth-like planet. See fig. [2]caption for description. 
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Figure 8. Orbital parameters evolution for the Keplerian models 
and the Jovian planet. See fig. [2] caption for description. 



viscosity regimes, the position of the shock front is generally 
closer to the central accreting object. 

Together with the above specified simulations, we per- 
formed a numerical test with our code, without the planet 
and with an initially totally Keplerian distribution of veloc- 
ities. The aim of this test is to analyse the survival time of 
the Keplerian kinematics before consistent radial perturba- 
tions comes out. Results are shown in figure IT71 where the 
ratio "radial speed/Keplerian speed" is plotted as a function 
of the radial distance to the central star. Four snapshots are 



Table 3. Linear best-fit parameters for the specific angular mo- 
mentum time evolution. For some cases, the linear fit is calculated 
over a limited time interval, during which the angular momentum 
evolution is linear. The linear fits are not present for the Keplerian 
models since there is no time interval during which the evolution 
of the angular momentum is linear. Migration times estimated 
with the formula l z /l z are reported in the last two columns: in 
the fourth column the dimensionless time is reported, while in 
the last column the same time is converted in years by choosing 
1 Mq and 2 Rq as physical reference units. 



provided at different times, the first at time zero, and the 
last at time — 3000 (corresponding to three Keplerian or- 
bits at r = 100). A spontaneous breaking of the initially 
Keplerian kinemetics is observed, due probably to the con- 
tribution of pressure forces and to the collisional nature of 
our model. 



4 DISCUSSION 

In a standard disc model, the relation connecting the disc 
thickness H, the disc radius r is: H/r = c s /vk ep i, where 
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Figure 10. Data and exponential best-fit for angular momentum 
in j = 18 models 



Figure 13. Data and exponential best-fit for angular momentum 
in Keplerian models 
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Figure 11. Data and exponential best-fit for angular momentum 
in j = 36 models 



c s is the sound velocity and Vk ep i is the Keplerian veloc- 
ity at the disc outer edge. Such relation correlates the disc 
geometry to a tight Keplerian kinematic condition. In our 
Keplerian disc model, the velocity of particles at the injec- 
tors positions is vtepi = 0.55 in the units specified in section 
12.11 c a — 0.05 being the sound speed, so that we obtain 
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Figure 12. Data and exponential best-fit for angular momentum 
in j = 54 models 
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Figure 14. Wave propagation in the j = 18 accretion disc. The 
four snapshots are ordered according to increasing time, start- 
ing from the bottom-left square in clock-wise order. The arrows 
indicate the position of the shock front. 



H/r ~ 0.05/0.55 ~ 0.09. Generally this evaluation applies 
to a thin disc, where the hypothes is of vertical hydrost atic 
equilibrium is taken into account (|Terquem et aLlboOOT ). If 
this relation is translated also to sub-Keplerian modelling, 
it would give H/r ~ 0.12, 0.18, 0.36 for j = 54, 36, 18 respec- 
tively. However, although these results could give some indi- 
cation, at least as a rough evaluation, these disc spreads are 
not strictly correct whenever turbulence and sub-Keplerian 
conditions are taken into account. The results of this work 
deal with protoplanet migration in a protoplanetary accre- 
tion disc around a pre-main sequence star, according to the 
hypothesis that the accretion process is not strictly Kep- 
lerian. The Keplerian constraint is often used in computa- 
tional models of accretion proto-planetary accretion discs on 
the basis of observational features that can be reproduced 
with Keplerian discs. It is sustainable if pressure and other 
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Figure 15. Wave propagation in the j = 36 accretion disc. The 
four snapshots are ordered according to increasing time, start- 
ing from the bottom-left square in clock-wise order. The arrows 
indicate the position of the shock front. 



Figure 16. Distribution of the gas particles around a Jupiter-like 
protoplanet in the initially Keplerian accretion disc. This snap- 
shot corresponds to a At ~ 3300 after the artificial formation of 
a purely Keplerian velocity distribution, which corresponds to a 
bit more than 3 orbits of the protoplanet. 



forces are not able to significantly modify an initially Kep- 
lerian distribution of velocities. We decided to explore the 
consequences of the presence of sub-Keplerian flows as a 
parametric study on the initial angular momentum, releas- 
ing the Keplerian constraint. 

Our models are technically considered as physically invis- 
cid, since the y are based on the Euler equations. Artifi- 
cial viscosity l|Monag han 1992]) is necessarily introduced in 
our models to resolve shocks numerically, to avoid spuri- 
ous heating and to handle fluid discontinuities. Artificial 
viscosity vanishes when the limit value of the particle in- 
terpolation domain goes to zero. The linear o>sph and 
quadratic f3s ph artificial viscosity terms (usually ~ 1 and 
sometimes, in som e specific cases, < 1) a re chosen t o be 1 
and 2, respect i vely. [Meglicki et all (ll993T l ; iDrimmelll (|l99rj ) ; 
iMurravl (|l996t ); IOkazaki et alj (|2002h demonstrated that the 



v = assCsH = c s h 



ass = h/H 



(3) 



linear component of the artificial visc osity, in t he co ntin- 
uum limit, yie l ds a v iscous shear force. IMurravl (|l996t ) and 
lOkazaki et "all l|2002l) explicitly formulated such an artificial 
viscosity contribution in momentum and energy equations, 
founding an analogy between the shear viscosity generated 
by the linear artific i al vis cosity term and the well-known 
Shak ura fc Sunvaevl (1 19751 ) shear viscosity in the contin- 
uum limit. In the SPH method, the quadratic (/3sph, von 
Neumann-Rychtmyer like viscosity) artificial viscosity term 
is used in order to handle structures. An analytical formula- 
tion, describi ng the numerical arti ficial viscosity coefficient, 
is reported in lMolteni et al.l (| 199lT ) : v num = c s h, where c s is 
the sound velocity. Therefore, in our modelling, any radial 
transport process is due to the role of the artificial viscosity. 
In order to estimate what ass value the numerical viscos- 
ity ref ers to, we can use the equality (see lLanzafame et all 
(2008) and references therein for a more detailed discussion 
of this point): 



Using the above estimated H/r values we obtain: ass 



2.5 x 10" 



1.9 x 10" 



1.3x10"^; 6.4 x 10"^ for Keple- 



rian and j = 54, 36, 18 respectively. 

In order to allow a comparison between the effects pro- 
duced by our implementation of artificial viscosity and the 
effects produced by physical shear viscosit y, we performed a 
test comparing a simulated spreading ring (|Speith fc Riffertl 
1 19991 : ISpeith fc Klevlliool) with its approximate analytical 
solution. As in the above specified papers we considered an 
initial ring of 40000 particles with a distribution of density 
given by the analytical solutio n obtained putting the di- 
mensionless time r = Y2vtjR\ (|Speith fc Kiev! [20031 ) equal 
to 0.018. The initial ring radius Rq is chosen to be equal to 
the stellar radius, while the kinematic viscosity is estimated 
as v ~ c 3 h, where h w as set equal to 0.09 as specified by 
ISpeith fc Riffertl i| 19991 ). 

Graphs showing the results of this test, at various times 
(specified by tau values) are shown in figure 1181 wh i ch ca n 
be directly compared to figure 1 in ISpeith fc Klevl {2003). 
Beware that this comparison, while allowing to show the 
sh ear effects produced by artificial viscosity, has its limits. 
In lSpeith fc Klevl (|2003) the disc is totally pressure-less and 
viscosity is treated via the viscous stress tensor. In our case, 
even though we artificially removed the gas pressure, we 
could not remove the pressure contribution due to the arti- 
ficial viscosity and we did not introduced a direct treatment 
of shear viscosity (our code is based on the Euler equations) . 
Moreover, the relation [3] links the artificial viscosity and the 
kinematic viscosity to the value of the speed of sound c s , so 
that the gas cannot be considered pressure- less. The main 
difference between o ur test results and those obtained by 
ISpeith fc Klevl (|2003l ) is the absence, in our case, of spiral 
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Figure 17. Distribution of the ratio "radial speed/Keplerian 
speed" as a function of radial distance for an initially totally 
Keplerian model. The four graphs show snapshots obtained at 
various times, considering that time = 1000 is one orbital period 
at a distance r = 100. 



structures. 

Concerning the difference between the faster migration in 
our sub-Keplerian models and slower migration in Keple- 
rian ones, we can notice that: in sub-Keplerian models the 
most important contribution to downward migration seems 
to come from momentum transfer between the gas particles 
and the protoplanet; the sub-Keplerian j x value induces an 
average decrease of planet l z , until the accretion disc cen- 
trifugal barrier is reached; an initially Keplerian distribu- 
tion of gas particles velocities, like in our Keplerian models, 
makes this effect less important, particularly with the Jovian 
planet. Models kepLa and kepljb show in fact a slow out- 
ward migration, despite the opening of a density gap. This 
result is different from classical views on the subject, which 
connect the opening of a density gap to a viscosity driven 
downward migration (type II migration), but without phys- 
ical viscosity (our models) the result can be different. 
Some people dPapaloizou et aL1l2006l ) suggest a type III mi- 
gration mechanism, where the main migration "engine" is 
provided by material flowing through the co-orbital region. 
This material would change orbit (and angular momentum) 
while passing trough the co-orbital, exerting a torque onto 
the protoplanet. this torque can be evaluated trough mo- 
mentum conservation. This effect might actually be present 
in our models, but it's not the main mechanism, at least in 
our sub-Keplerian models. In our case, momentum conserva- 
tion is directly applied to correct the protoplanet kinematics 
when disc particles are captured by the planet. A rough es- 
timate of the specific angular momentum change due to the 
capture of one disc particle can be done as follows: let I, j, 
M and m be the specific angular momentum and mass of 
the protoplanet and of the gas particle respectively; angu- 
lar momentum conservation after particle capture leads to 




Figure 18. Evolution of a ring of partic les with an initial dis - 
tribution of 40000 particles as specified in ISpeith fe Kiev! tOO^I . 
X and Y positions are in units of Ro and snapshots at various r 
values are shown (see text). 



Al = m/(M + m)(j — I) ~ m/M(j - I) for m/M « 1; in 
our case (planet at an initial distance of 100 in a circular Ke- 
plerian orbit) Unitiai = 62.8; referring to our 36c model we 
have m/M — 10 -4 and j = 36, so that after a single capture 
Al ~ -2.7 x 10~ 3 . About 10 s particles are captured during 
our simulation in model 36c before the angular momentum 
is reduced to about one half of the initial value (not all of 
them will have preserved the j = 36 before capture), so that 
this capture mechanism is perfectly able to account for our 
rapid migration. 

In our Keplerian models particle captures can both increase 
or decrease the protoplanet angular momentum, and a flow 
of particles through the co-orbital region is visible (figure 
I16p . so type III mechanism might play a role. Notice how- 
ever that the above specified flowing particles are moving 
towards smaller orbits, so they are mostly loosing angular 
momentum. This might actually be the reason for our slowly 
increasing angular momentum of the protoplanet in fcepLa, b 
models. But, being particle captures present in all our mod- 
els, we cannot have pure type III migration anyway. 

The results also lead our attention to the secondary role 
of the presence of a "pseudo-atmosphere" surrounding the 
protoplanets with the aim to discover whether such cumu- 
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lative role could be also relevant. 

The impact of the eventual presence of a pseudo-atmosphere 
in the protoplanets is visible in some models, particularly 
those with the less massive protoplanets (look at the differ- 
ence between models 54c and 54d), but never appears crucial 
in determining the main results. The slightly faster migra- 
tion rate shown by models without the pseudo-atmosphere 
was actually expected, because the pseudo-atmosphere adds 
a repulsive pressure interaction between the protoplanet and 
the gas, that should hinder gas particle capture by the pro- 
toplanet. Since part of the migration effect in our models 
is due to momentum transfer between the infalling particles 
and the protoplanets when the former are captured, a lower 
capture rate might bring a lower migration rate. 
Particularly interesting is the result reported by the c, d 
models, where the radial migration breaks down close to 
~ 15; 20; 45 stellar radii, for j — 18; 36; 54 respectively. We 
believe that this effect is mainly due to the interaction of 
the protoplanet with the propagation of outward gas waves 
coming from inner zones (figures [T4"lll5p . where the centrifu- 
gal barrier is located in our sub-Keplerian models. These 
results show that, despite the sub-Keplerian and more real- 
istic conditions, planetary migration could even halt. Sub- 
sequent disc extinction could freeze the protoplanet orbital 
parameters. Of course, due to a much higher inertia, Jovian 
planets evolve much more slowly. 

An analytical expression giving us an idea of the protoplanet 
evolution, as a function of the model initial conditions, is ex- 
pressed by the specific angular momentum fit equations I H2I 
and by their temporal first derivatives. 

Even if it is not the main topic of this paper, let's draw 
some comments on eccentricity behaviour shown in the bot- 
tom rectangles of figures I2"l9l For Jupiter-like planets the or- 
bits remain almost circular, with values around 10 -5 , while 
for Earth-like planets some eccentricity enhancing effect is 
visible in sub-Keplerian models, but the behaviour is quite 
complex and a more stable value is obtained in late migra- 
tion phases, with values in the 0.01 — 0.02 range. 
For a compar ison with other resu lts we refer to the study 
performed bv lSchafer et alj |2004) since they used SPH to 
simulate accretion disc hydrodynamics, even if some differ- 
ences have to be remarked: first, they use a different nu- 
merical method to avoid particle mutual penetration, since 
they rely on XSPH with a bulk viscosity contribution to 
the viscous stress tensor; second, they simulate ph ysical 
viscous discs according to IShakura fc Sunvaevl (1 19731) pre- 
scription whilst our discs are inviscid; third, they use an 
isothermal equation of state whilst we rely on a classical 
ideal gas equation with variable thermal energy per unit 
mass; fourth, their initial velocity distribution is Keplerian. 
They estimate a decrease of the planet distance from the 
the star of about 1 AU (from 5.2 AU to 4.1 AU) in about 
30000 years for a Jovian planet, which brings to a rate of 
Da/Dt ~ 3.3 x 10~ 5 AU/y. They also conclude that migra- 
tion is braked by disc mass loss, which means a lowering of 
the disc density and of the gravitational interaction. 
Inserting in our dimensionless model the following reference 
physical units: central star mass Mp = IM q; stellar radius 
from typical T-Tauri stars ()Bertoudll989h R ~ 27?©; the 
basic unit for time is To ~ 2.8 x 10 4 s. Our sub-Keplerian 
models with j = 18; 36; 54 give then migration rates of 
Da/Dt ~ 8 x 10~ 5 ;5.3x~ 5 ;3.4 x 10 -5 AU/y respectively. 



We notice that our j — 54 (the less sub-Keplerian model) 
migration rate is similar to th e average migration rate ob- 
tained bv lSchafer et""ai1 [|2004h . even if our starting distance 
is about 1 AU with the above specified choices for physical 
units. More important, the braking effect is due, in our simu- 
lations, to a completely different mechanism based on turbu- 
lence, outward propagation of bouncing waves and changing 
of the average momentum transferred to the planet during 
gas particle captures, and we do not observe a significant 
lowering of the disc density. 
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